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Abstract 


The  evolution  of  laminar  vortices  was  calculated  using  the  Reynolds-averaged  Navier- 
Stokes  solver  ANSYS  CFX.  Vortices  of  different  strengths  were  modelled  on  unstruc¬ 
tured  grids  using  different  element  types  and  different  node  densities,  then  compared 
to  an  approximate  solution  of  the  Navier-Stokes  equation.  Grids  using  hexahedral 
or  prismatic  elements  aligned  with  the  vortex  yielded  significantly  more  accurate 
solutions  than  tetrahedral  grids,  especially  when  the  circulating  velocity  was  much 
smaller  than  the  convection  velocity.  Over-diffusion  of  the  circumferential  velocity  is 
the  major  source  of  error  in  these  solutions. 

ANSYS  CFX  was  also  used  to  determine  the  effect  of  grid  density  on  a  turbulent 
vortex,  though  in  this  case  there  is  no  analytic  solution  to  compare  with.  As  with  the 
laminar  vortex,  grids  using  hexahedral  or  prismatic  elements  aligned  with  the  vortex 
converge  much  more  rapidly  with  cell  size  than  tetrahedral  grids,  suggesting  that 
they  are  more  accurate.  In  these  solutions  the  over-diffusion  of  turbulent  viscosity  by 
artificial  viscosity  is  an  important  source  of  error. 

The  solutions  are  used  to  make  recommendations  for  the  density  of  unstructured  grids 
required  for  accurate  predictions  of  the  pressure  in  the  core  as  the  vortex  evolves  over 
hundreds  of  core  diameters. 


Resume 


On  a  calculc  revolution  de  tourbillons  laminaires  en  utilisant  le  logicicl  de  resolution 
de  1’ equation  de  Navier-Stokes  a  moyenne  de  Reynolds  ANSYS  CFX.  On  a  modelise 
des  tourbillons  d’intensites  differentes  sur  des  maillages  non  structures  en  utilisant 
different  s  types  d’clcments  et  differentes  densites  nodales,  puis  en  procedant  a  une 
comparaison  avec  une  solution  approximative  de  l’equation  de  Navier-Stokes.  Les 
maillages  a  elements  hexahedraux  prismatiques  alignes  avec  lc  tourbillon  ont  donne 
des  solutions  nettement  plus  exactes  que  les  maillages  tetrahedraux,  particulicrement 
quand  la  vitesse  de  circulation  etait  beaucoup  plus  petite  que  la  vitesse  de  convection. 
La  surdiffusion  de  la  vitesse  circonferenticlle  est  la  principale  source  d’erreur  dans  ces 
solutions. 

On  a  egalement  utilise  ANSYS  CFX  pour  determiner  l’effet  de  la  densite  du  maillage 
sur  un  tourbillon  turbulent,  bien  qu’il  n’y  ait  pas  dans  ce  cas  de  solution  analy- 
tique  permettant  de  faire  une  comparaison.  Comme  avec  lc  tourbillon  laminaire,  les 
maillages  a  elements  hexahedraux  ou  prismatiques  alignes  avec  le  tourbillon  convergent 
beaucoup  plus  rapidement  avec  la  taille  des  mailles  que  les  maillages  tetrahedraux, 
ce  qui  porte  a  croire  qu’ils  sont  plus  precis.  Dans  ces  solutions,  la  surdiffusion  de  la 
viscosite  turbulcnte  par  la  viscosite  artificielle  est  une  importante  source  d’erreur. 
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Les  solutions  sont  utilisees  pour  faire  des  recommandations  concernant  la  densite 
des  maillages  non  structures  qui  est  requise  pour  obtenir  des  previsions  exactes  de 
la  pression  au  coeur  du  tourbillon  quand  celui-ci  couvre  des  distances  de  1’ordre  de 
plusieurs  centaines  de  diametres  du  coeur. 
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Executive  summary 

RANS  calculations  of  the  evolution  of  vortices  on 
unstructured  grids 

David  Hally;  DRDC  Atlantic  TM  2009-052;  Defence  R&D  Canada  -  Atlantic; 
November  2009. 

Background:  The  Maritime  Asset  Protection  Section  at  DRDC  Atlantic  is  applying 
Reynolds-averaged  Navier-Stokes  (RANS)  flow  solvers  to  predict  the  flow  around  ship 
propellers  with  the  goal  of  being  able  to  model  the  noise  generated  from  propeller 
cavitation.  Of  particular  interest  is  the  cavitation  in  the  vortex  produced  at  the 
propeller  tip.  This  form  of  cavitation  is  often  the  first  to  appear.  It  cannot  be 
predicted  by  simpler  flow  solvers  (panel  methods)  available  at  DRDC  Atlantic. 

In  order  to  predict  the  cavitation,  an  accurate  estimate  of  the  pressure  in  the  vortex 
core  is  needed.  One  question  of  interest  is  how  dense  the  computational  grid  must 
be  in  the  vortex  core  in  order  to  model  the  pressure  adequately.  The  answer  will 
have  an  important  effect  on  the  gridding  strategy  used  when  preparing  a  propeller 
flow  calculation.  This  report  examines  the  effects  of  grid  type  and  density  on  the 
predictions  of  the  pressures  in  the  cores  of  long  rectilinear  vortices. 

Results:  The  evolution  of  laminar  and  turbulent  vortices  was  calculated  using  the 
RANS  solver  ANSYS  CFX  from  ANSYS,  Inc.  Vortices  of  different  strengths  were 
modelled  on  unstructured  grids  using  different  element  types  and  different  node  den¬ 
sities.  The  solutions  for  laminar  vortices  were  compared  to  an  approximate  solution 
of  the  Navier-Stokes  equation. 

It  was  found  that  significantly  more  accurate  solutions  were  obtained  when  the  grid 
cells  were  aligned  with  the  direction  of  the  vortex,  especially  when  the  circulating 
velocity  was  much  smaller  than  the  convection  velocity.  In  particular,  grids  using 
only  tetrahedra  (in  effect,  no  alignment  with  the  flow  at  all)  cause  the  vortex  to 
dissipate  rapidly  unless  very  large  numbers  of  cells  are  used. 

Significance:  These  calculations  indicate  that  although  ANSYS  CFX  can  predict 
vortex  flows  accurately  using  unstructured  grids,  the  size  of  the  grids  may  become 
very  large  for  the  calculation  of  propeller  tip  vortices.  The  estimates  of  minimum  grid 
density  will  be  used  when  formulating  gridding  strategies  for  vortices  on  propellers. 

Future  work:  The  calculations  reported  will  be  used  to  develop  a  strategy  for  gen¬ 
erating  grids  capable  of  resolving  propeller  tip  vortices  adequately. 
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RANS  calculations  of  the  evolution  of  vortices  on 
unstructured  grids 

David  Hally ;  DRDC  Atlantic  TM  2009-052 ;  R  &  D  pour  la  defense  Canada  - 
Atlantique ;  novembre  2009. 

Contexte  :  La  Section  de  la  protection  des  biens  maritimes  de  RDDC  Atlantique 
utilise  des  logiciels  de  resolution  de  l’equation  de  Navier-Stokes  a  moyenne  de  Rey¬ 
nolds  (RANS)  pour  prevoir  l’ecoulement  autour  des  helices  de  navire  dans  le  but  de 
modeliser  le  bruit  de  la  cavitation  cause  par  les  helices.  La  cavitation  dans  le  tour- 
billon  cree  a  l’extremite  des  helices  presente  un  interet  particulier.  Cette  forme  de 
cavitation  est  souvent  la  premiere  a  apparaitre.  Les  logiciels  de  calcul  d’ecoulement 
(par  la  methode  de  la  plaque  pleine)  presentement  disponibles  a  RDDC  Atlantique 
ne  permettent  pas  de  la  prevoir. 

Une  estimation  exacte  de  la  pression  au  coeur  du  tourbillon  est  necessaire  pour  pre¬ 
voir  cette  cavitation.  L’une  des  questions  a  examiner  est  la  densite  que  doit  avoir 
le  maillage  de  calcul  dans  le  coeur  du  tourbillon  pour  que  la  pression  soit  modelisee 
de  fagon  adequate.  La  reponse  aura  une  repercussion  importante  sur  la  strategic  de 
maillage  utilisee  dans  le  calcul  de  l’ecoulement  autour  de  l’helice.  Le  present  rapport 
examine  les  effets  du  type  et  de  la  densite  de  maillage  sur  les  previsions  de  la  pression 
au  coeur  de  tourbillons  rectilignes  longs. 

Resultats  :  On  a  calcule  revolution  de  tourbillons  laminaires  et  turbulents  en  uti- 
lisant  le  logiciel  de  resolution  de  l’equation  de  Navier-Stokes  a  moyenne  de  Reynolds 
ANSYS  CFX  de  la  societe  ANSYS  Inc.  On  a  modelise  des  tourbillons  d’intensites 
differentes  sur  des  maillages  non  structures  en  utilisant  differents  types  d’elements 
et  differentes  densites  nodales.  Les  solutions  pour  les  tourbillons  laminaires  ont  ete 
companies  a  une  solution  approximative  de  l’equation  de  Navier-Stokes. 

On  a  obtenu  des  solutions  nettement  plus  exactes  quand  les  mailles  etaient  alignees 
avec  la  direction  du  tourbillon,  particulierement  quand  la  vitesse  de  circulation  etait 
beaucoup  plus  petite  que  la  vitesse  de  convection.  En  particulier,  a  moins  d’utiliser 
un  grand  nombre  de  mailles,  les  tourbillons  subissent  une  dissipation  rapide  dans 
les  maillages  qui  n’utilisent  que  des  tetrahedres  (absence  totale  d’alignement  avec 
l’ecoulement). 

Importance  :  Ces  calculs  indiquent  que,  bien  que  le  logiciel  ANSYS  CFX  puisse 
prevoir  les  ecoulcments  de  tourbillons  avec  exactitude  au  moyen  de  maillages  non 
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structures,  1c  nombre  de  mailles  necessaire  peut  devenir  tres  grand  pour  le  calcul  des 
tourbillons  a  l’extremite  des  helices.  On  utiliscra  les  estimations  des  densites  mini¬ 
males  des  maillages  dans  la  formation  des  strategies  de  maillage  pour  les  tourbillons 
causes  par  les  helices. 

Perspectives  :  Les  calculs  publies  seront  utilises  pour  claborer  une  strategic  de 
generation  de  maillages  permettant  de  calculcr  adequatement  les  tourbillons  aux  ex- 
tremites  des  helices. 
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1  Introduction 


The  Maritime  Asset  Protection  and  the  Warship  Performance  Sections  at  DRDC 
Atlantic  apply  Reynolds-averaged  Navier-Stokes  (RANS)  solvers  to  the  prediction 
of  vortices  generated  by  ships,  submarines  and  propellers.  Important  applications 
are  the  prediction  of  tip  and  leading  edge  vortex  cavitation  on  propellers,  of  keel 
and  body  vortices  which  can  influence  propeller  efficiency  and  underwater  vehicle 
maneuvering  characteristics,  and  of  low  aspect  ratio  appendage  wake  vortices  which 
can  compromise  the  performance  of  downstream  control  surfaces  (for  example,  the 
influence  of  a  submarine  sail  on  its  tailplanes). 

As  a  preliminary  test,  Hally  and  Watt[l]  calculated  the  evolution  of  a  laminar  vortex 
on  structured  grids  using  the  RANS  solvers  TRANSOM[2],  developed  at  DRDC  At¬ 
lantic,  and  CFX-TASCflow[3],  developed  by  AEA  Technology  Ltd  (now  ANSYS,  Inc.). 
There  is  an  approximate  solution  to  the  Navier-Stokes  equations  which  describes  the 
vortex  making  it  a  good  verification  case  for  the  flow  solvers.  This  document  describes 
similar  calculations  on  unstructured  grids  using  the  RANS  solver  ANSYS  CFX. 

Four  different  types  of  grid  were  constructed  using  different  combinations  of  hexahe- 
dral,  prism  and  tetrahedral  elements.  Of  interest  is  the  effect  of  the  grid  density  on 
the  solution:  for  each  type  of  grid,  how  many  cells  are  necessary  in  the  vortex  core  to 
model  the  pressure  adequately?  The  dependence  on  the  grid  will  have  an  important 
effect  on  the  gridding  strategy  used  when  modelling  propeller  vortices.  The  depen¬ 
dence  on  grid  density  was  tested  using  the  laminar  vortex  as  well  as  predictions  of  a 
turbulent  vortex  for  which  there  is  no  analytical  solution. 


2  The  evolution  of  a  laminar  vortex 


2.1  Theory 


In  a  convection  dominated  flow,  in  a  cylindrical  coordinate  system  (see  Fig.  1),  an 
approximate  steady  solution  to  the  Navier-Stokes  equations  is: 


vx  =  U;  vr  =  0;  ve 


C 

2tt  r 


P  = 


pC2U 
327 t2ux 


/(*) 


(1) 


where 


r2U 

Avx' 


/(*) 


3-2/ 


2 

dy 


(2) 


This  is  a  vortex  with  centre  on  the  line  r  =  0,  circulation  C,  fluid  density  p  and 
kinematic  viscosity  v.  The  vortex  decays  as  x  increases.  Annex  C  shows  that  the 
solution  is  a  good  approximation  provided  that 


v  -C  0.2  rcll;  (3  <9 


(3) 
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Figure  1:  The  coordinate  system 


where  rc  is  the  radius  of  the  vortex  core,  defined  to  be  the  location  of  maximum  Vg\ 


rc 


(4) 


and 


P  = 


C 

UTr, 


(5) 


with  rCi  equal  to  rc  in  the  upstream  plane.  The  parameter  (3  is  a  non-dimensional 
measure  of  the  strength  of  the  vortex. 


The  maximum  value  of  Vg  is 


C 

^0 max  ~  0.114 

TV 


(6) 


More  rehned  approximations  for  the  axial  and  radial  velocities,  vx  and  vr,  are  given 
in  Annex  C. 


Eq.  (4)  implies  that  the  core  radius  of  rci  will  double  to  2 rci  in  a  distance 


r  1.6Ur2r 
v 


(7) 


Over  the  same  distance  the  maximum  velocity  drops  by  half  and  the  pressure  drop 
in  the  core  decreases  by  one  quarter. 


2.2  RANS  calculations 

ANSYS  CFX  was  used  to  calculate  the  evolution  of  a  laminar  vortex  in  water  at 
standard  temperature  and  pressure.  In  all  cases  the  grid  was  bounded  by  a  cylinder  of 
radius  R  =  20  metres  with  axis  extending  from  xq  =  1000  metres  to  X\  =  2000  metres. 
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The  convection  velocity  was  set  by 


v 


U 


10-3m 


(8) 


and  the  circulation  was  set  to  C  =  U  x  0.1m,  U  x  lm  or  U  x  10m.  The  lowest  value 
corresponds  to  the  earlier  calculations  by  Hally  and  Watt.  The  three  values  of  C  used 
in  the  calculations  correspond  to  /3  values  of  0.0446,  0.446  and  4.461. 


The  upstream  plane  was  an  inlet  with  velocity  set  using  Eq.  (1).  The  outer  boundary 
was  a  free  slip  wall.  The  downstream  plane  was  an  outlet  with  average  normal  velocity 
set  to  U. 


The  ANSYS  CFX  “High  Resolution  Scheme”  for  discretization  of  the  convection  terms 
was  used  because  it  was  felt  that  this  would  be  the  most  practical  scheme  for  the  more 
complex  flow  around  a  propeller.  The  effect  of  this  scheme  relative  to  the  second  order 
accurate  scheme  is  discussed  in  Sec.  2.8. 


Four  different  types  of  grid  were  used: 

1.  a  structured  grid  consisting  of  hexahedra  except  for  those  elements  along  the 
centreline  which  were  triangular  prisms; 

2.  a  hybrid  grid  with  an  inner  portion  using  hexahedra  and  an  outer  portion  using 
triangular  prisms  aligned  with  the  flow; 

3.  a  grid  consisting  of  triangular  prisms  aligned  with  the  flow;  and 

4.  a  grid  consisting  of  tetrahedra. 

For  each  grid  type,  calculations  were  made  using  five  grids  of  differing  node  density 
near  the  vortex  core.  Each  grid  is  described  in  detail  in  the  following  sections. 

All  the  grids  were  generated  using  the  scripting  language  Pointwise  Glyph  and  the 
program  Pointwise [4]  from  Pointwise,  Inc. 


2.3  The  hexahedral  grids 


Let  Nc  be  the  number  of  nodes  across  the  core  of  the  vortex.  The  hexahedral  grids 
had  nodes  of  the  form  {(ri,9j,Xf.)  ■  i  G  [1  ,Nr\;j  €  [1,  Ng];k  e  [1,  Nx]}  with 


_  27t(z  —  1)  _  (k-  l)(x1-x0) 

i  -  -  1  X  b  -  Xf)  \  - 

3  Ne-l  ’  Nx  —  1 


(9) 


In  all  cases  Ng  =  64  and  Nz  =  101.  The  r*  were  calculated  using  a  tanh  distribution 
(see  Hally[5])  with  the  spacing  at  r  =  0  set  to  s0  =  2 rc/Nc  and  the  spacing  at  r  =  R 
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set  to  sr  =  R/ 25.  The  value  of  Nr  was  set  so  that  the  tanh  distribution  approximated 
a  geometric  distribution  as  closely  as  possible:  i.e.  the  ratio  of  the  cell  sizes  in  the 
radial  direction  increased,  as  nearly  as  possible,  in  a  geometric  progression: 


Nr 


1  R\n(sR/so) 

Sr  —  So 


(10) 


where  [x]  denotes  the  closest  integer  to  x. 


The  left  side  of  Fig.  2  shows  the  grid  on  the  upstream  face  of  the  cylinder  for  the  case 
with  Nc  =  10. 


2.4  The  hybrid  grids 

The  hybrid  grids  were  generated  by  first  making  a  hybrid  grid  on  the  upstream  face 
of  the  cylinder.  This  grid  was  then  extruded  downstream  using  Nx  translation  steps 
parallel  to  the  cylinder  axis.  With  Nx  =  101,  this  resulted  in  prismatic  elements  each 
of  length  10  metres. 

The  hybrid  grid  on  the  upstream  face  was  generated  in  three  steps. 

1.  First  a  structured  block  two  nodes  deep  was  generated  around  the  circumference 
of  the  cylinder.  The  purpose  of  this  block  was  to  improve  the  application  of  the 
CFX  boundary  conditions  at  the  outer  boundary. 

2.  Then  a  square  structured  block  with  sides  of  length  6 rc  was  created  to  cover 
the  core  of  the  vortex.  The  node  separation  in  this  block  is  rc/Nc. 

3.  Finally  an  unstructured  block  was  created  in  the  region  between  the  previous 
two  blocks.  The  node  separation  increased  gradually  from  the  inner  to  the  outer 
boundary. 

The  right  side  of  Fig.  2  shows  the  grid  on  the  upstream  face  of  the  cylinder  for  the 
case  with  Nc  =  10. 


2.5  The  grids  with  prismatic  elements 

The  grids  with  prismatic  elements  were  generated  by  first  making  an  unstructured 
grid  on  the  upstream  face  of  the  cylinder.  This  grid  was  then  extruded  downstream 
in  the  same  manner  as  the  hybrid  grids. 

The  unstructured  grid  on  the  upstream  face  was  generated  in  a  manner  similar  to  the 
hybrid  grids  except  that  the  inner  block  was  an  unstructured  block  of  radius  rc  with 
mean  node  separation  rc/Nc. 

Fig.  3  shows  the  grid  on  the  upstream  face  of  the  cylinder  for  the  case  with  Nc  =  10. 
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Figure  2:  The  upstream  face  of  the  hexahedral  grid  (left)  and  hybrid  grid  (right) 
with  Nc  =  10.  The  red  circle  marks  the  size  of  the  vortex  core. 


Figure  3:  The  upstream  face  of  the  prismatic  and  tetrahedral 
grids  with  Nc  =  10.  The  red  circle  marks  the  size  of  the  vortex 
core. 
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2.6  The  grids  with  tetrahedral  elements 

The  grids  with  tetrahedral  elements  were  generated  in  four  steps. 

1.  First  a  structured  block  two  nodes  deep  was  generated  on  the  surface  of  the 
cylinder.  As  with  the  hybrid  and  prismatic  grids,  the  purpose  of  this  block 
was  to  improve  the  application  of  the  CFX  boundary  conditions  at  the  outer 
boundary. 

2.  An  inner  cylinder  of  radius  rc  was  generated  and  its  surface  gridded  with  a  mesh 
of  triangles  with  mean  edge  length  of  rc/Nc. 

3.  The  interior  of  the  inner  cylinder  was  gridded  using  tetrahedra  with  mean  edge 
length  of  rc/Nc. 

4.  The  region  between  the  inner  cylinder  and  the  outer  structured  block  was  grid¬ 
ded  using  tetrahedra  whose  edge  lengths  increased  gradually  with  radius. 

The  grid  on  the  upstream  face  of  the  tetrahedral  grids  is  similar  to  that  of  the  pris¬ 
matic  grids:  see  Fig.  3. 

It  should  be  noted  that  the  tetrahedral  grids  have  far  more  elements  than  the  hexa- 
liedral,  hybrid  or  prismatic  grids.  This  is  because  the  tetrahedra  have  edges  that  are 
roughly  uniform  in  size  in  all  directions,  while  in  the  other  two  grids  the  edges  aligned 
with  the  cylinder  axis  are  much  longer  than  the  other  edges.  Table  1  gives  estimates 
of  the  number  of  nodes  and  elements  for  each  of  the  grids  as  well  as  the  number 
of  nodes  and  elements  in  the  vortex  core.  The  latter  two  number  increase  roughly 
linearly  with  Nc  for  the  hexahedral  grids,  quadratically  with  Nc  for  the  hybrid  and 
prismatic  grids,  and  as  for  the  tetrahedral  grids. 


6 


DRDC  Atlantic  TM  2009-052 


Table  1:  Estimates  of  the  number  of  nodes  and  elements  in  each  of  the  grids. 


Nc 

No.  Nodes 
total 

No.  Elements 
total 

No.  Nodes 
in  core 

No.  Elements 
in  core 

Hexahedral 

3 

114,635 

113,400 

12,827 

12,600 

6 

165,539 

163,800 

19,190 

18,900 

10 

210,080 

207,900 

31,916 

31,500 

16 

254,621 

252,000 

51,005 

50,400 

22 

292,799 

289,800 

70,094 

69,300 

Hybrid 

3 

52,823 

83,900 

404 

700 

6 

122,513 

197,600 

2,727 

2,800 

10 

278,861 

449,600 

7,575 

7,900 

16 

644,986 

1,034,200 

19,695 

20,100 

22 

1,172,913 

1,874,400 

37,875 

38,000 

Prismatic 

3 

50,601 

87,600 

1,919 

2,600 

6 

80,497 

146,800 

4,343 

6,600 

10 

159,681 

303,600 

13,029 

22,400 

16 

316,635 

614,400 

24,947 

44,200 

22 

561,257 

1,098,800 

53,530 

98,800 

Tetrahedral 

3 

303,508 

1,548,449 

10,080 

37,327 

6 

484,929 

2,645,971 

51,014 

232,517 

10 

1,007,983 

5,803,958 

300,023 

1,551,813 

16 

2,057,450 

12,136,019 

948,935 

5,134,227 

22 

4,270,040 

25,525,534 

2,484,084 

13,841,015 
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2.7  Results  of  calculations 


It  is  important  to  note  that  the  CFX  boundary  conditions  are  incapable  of  capturing 
the  true  conditions  exactly  at  the  upstream  and  downstream  boundaries.  Therefore, 
very  close  to  the  upstream  and  downstream  planes  the  solutions  exhibit  higher  errors 
than  in  the  interior  of  the  cylinder.  To  avoid  the  boundary  effects  we  compare  in 
a  cross-stream  plane  at  x  —  xmid  =  1500  metres,  halfway  along  the  cylinder.  The 
results  are  also  compared  along  the  centreline  of  the  vortex,  r  =  0. 

Examples  of  the  results  for  the  prismatic  grid  with  [3  =  0.0446  are  shown  in  Figs.  4-6. 
Figs.  4  and  5  show  the  predicted  values  of  pressure  and  circumferential  velocity,  Vg, 
axial  velocity  deficit,  U  —  vXl  and  radial  velocity,  vr,  as  a  function  of  r/rCj  in  the 
plane  x  =  Xmid ■  Fig.  6  shows  the  pressure  and  axial  velocity  deficit  as  a  function  of 
{x  —  xq )/rCi  (i.e.  the  number  of  core  radii  downstream  of  the  initial  plane)  along  the 
centreline  (r  =  0).  Similar  plots  for  all  grid  types  and  values  of  (3  can  be  found  in 
Annex  A. 

To  generate  Figs.  4  and  5,  the  results  hie  from  the  calculation  was  loaded  into  the 
CFX  post-processor,  then  interpolated  to  points  in  the  plane  x  =  xmid,  then  written 
out  to  a  data  hie.  Each  datum  is  shown  by  plotting  its  value  (pressure  or  velocity 
component)  versus  the  distance  of  the  point  from  the  axis.  For  the  hexahedral  grid 
(see  Fig.  A.l,  for  example),  which  exhibits  complete  axisymmetry,  all  the  data  for  a 
single  calculation  should  collapse  to  a  line.  However,  for  the  unstructured  grids  in 
which  the  symmetry  is  broken,  the  data  exhibit  some  spread. 

In  the  figures,  the  solid  black  line  shows  the  values  predicted  using  Eq.  (1)  or  by 
the  refinements  described  in  Annex  C  and  given  by  Eqs.  (C.33)  and  (C.60).  When 
/ 3  =  0.0446  this  provides  an  accurate  estimate  of  the  true  solution;  all  numerical 
calculations  should  converge  to  the  black  line.  When  (3  =  0.446  (Figs.  A.4-A.6), 
the  approximations  used  in  Eq.  (1)  are  at  the  limits  of  their  validity;  all  numerical 
calculations  should  converge  close  to  the  black  line  but  small  discrepancies  should  be 
expected  within  the  vortex  core.  Also,  as  described  above,  for  large  r  the  axial  velocity 
deficit  U  —  vx  is  expected  to  be  slightly  negative  due  to  the  imposition  of  zero  radial 
velocity  on  the  outer  boundary.  When  (3  =  4.461  (Figs.  A.7-A.9),  the  approximations 
used  in  Eq.  (1)  are  no  longer  valid;  one  should  not  expect  the  numerical  predictions 
to  converge  to  the  black  line;  however,  it  is  retained  in  the  plots  to  make  it  easier  to 
compare  between  the  results  at  different  circulations. 

Because  none  of  the  boundary  conditions  used  in  the  calculations  (see  Sec.  2.2)  set 
a  value  for  the  pressure,  the  pressure  is  only  defined  relative  to  an  unknown  additive 
constant.  To  remove  the  ambiguity,  the  computed  pressure  values  were  offset  so  that 
the  pressure  at  (x,  r,  6)  =  ( xmid ,  R}  0)  was  equal  to  the  theoretical  value  predicted  by 
Eq.  (1). 
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n/(xA-n) 


Prismatic  Grid:  p  =  0.0446 


Prismatic  Grid:  p  =  0.0446 


Figure  4:  Laminar  vortex:  the  pressure  coefficient  (left)  and  circum¬ 
ferential  velocity  (right)  as  a  function  of  r/rCi  in  the  plane  x  =  xmid 
with  (3  =  0.0446. 


Prismatic  Grid:  p  =  0.0446  Prismatic  Grid:  p  =  0.0446 


0.0002 
0.00015 
0.0001 
5e-05 
0 

-5e-05 
-0.0001 
-0.00015 
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22  nodes  » 
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* 

Figure  5:  Laminar  vortex:  the  axial  velocity  deficit  (left )  and  radial 
velocity  (right)  as  a  function  of  r /rCi  in  the  plane  x  =  xmid  with 
(3  =  0.0446. 


Prismatic  Grid:  p  =  0.0446 


Prismatic  Grid:  p  =  0.0446 


(x-XoVci 


(x-x0)/rci 


Figure  6:  Laminar  vortex:  the  pressure  coefficient  (left)  and  axial 
velocity  deficit  (right)  as  a  function  of  ( x  —  x0)/rci  on  the  centreline 
with  f3  =  0.0446. 
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|5  =  0.0446  p  =  0.0446 


Figure  7:  The  pressure  (left)  and  circumferential  velocity  vg  (right) 
as  a  function  of  r  in  the  plane  x  =  xmui- 


Like  the  upstream  and  downstream  boundaries,  the  free  slip  wall  boundary  condition 
used  on  the  outer  boundary  is  also  not  strictly  accurate.  In  the  theoretical  solution 
the  radial  velocity  at  the  outer  boundary  is  small  but  positive.  Therefore  there  is  a 
small  loss  of  fluid  mass  through  the  outer  boundary  which  exactly  compensates  for  the 
difference  in  fluid  mass  convected  through  the  upstream  and  downstream  boundaries 
due  to  the  increase  in  axial  velocity  deficit.  The  axial  velocity  deficit  in  the  theoretical 
solution  rapidly  approaches  zero  for  large  r.  However,  the  free  slip  boundary  condition 
requires  that  there  be  no  flux  of  fluid  through  the  outer  boundary;  therefore  there  is 
no  loss  of  mass  through  the  outer  boundary  and  the  mass  convected  through  any  cross 
stream  plane  must  be  constant.  Therefore  the  integral  of  the  axial  velocity  deficit 
over  any  cross  stream  plane  must  be  zero.  This  causes  a  small  negative  shift  of  the 
calculated  axial  velocity;  at  large  r  the  calculated  axial  velocity  deficit  approaches  a 
small  negative  value,  not  zero. 

The  flow  with  (3  =  0.0446  corresponds  to  the  tests  run  earlier  by  Hally  and  Watt.  The 
current  results  for  this  flow  are  shown  in  Figs.  A.1-A.3.  For  this  case  the  hexahedral, 
hybrid  and  prismatic  grids  converge  very  well  even  when  the  number  of  nodes  across 
the  core  is  small.  These  results  conform  very  closely  to  those  reported  for  TRANSOM 
and  CFX-TASCflow  by  Hally  and  Watt.  However,  when  the  tetrahedral  grid  is  used 
the  results  are  very  poor;  even  with  22  nodes  across  the  core  the  solution  has  not  yet 
converged  adequately.  Fig.  7  shows  the  pressure  and  circumferential  velocity  predicted 
using  the  hexahedral,  hybrid  and  prismatic  grids  with  Nc  =  6  and  the  tetrahedral 
grid  with  Nc  =  22.  The  three  former  grids  give  much  more  accurate  results  despite 
the  significantly  smaller  numbers  of  nodes. 

The  primary  difference  between  the  tetrahedral  grid  and  the  other  three  grid  types  is 
that  the  convection  velocity,  vx,  which  is  the  dominant  portion  of  the  velocity  field, 
is  aligned  with  the  faces  of  the  prisms  and  hexahedra  (i.e.  it  is  either  normal  to  or 
parallel  to  a  face)  while  it  is  not  aligned  with  the  faces  of  the  tetrahedra.  This  suggests 
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P  =  4.46 


P  =  4.46 


Figure  8:  The  pressure  (left)  and  circumferential  velocity  vg  (right) 
as  a  function  of  r  in  the  plane  x  =  xmid  when  (3  =  4.461. 

that  artificial  viscosity  is  being  introduced  when  the  velocity  crosses  a  face  obliquely. 
Because  \vx\  \vg\,  the  artificial  viscosity  induced  by  vx  significantly  increases  the 
effective  viscosity  causing  vg  and  the  pressure  to  be  over-damped. 

Notice  that  with  the  hexahedral  grid  the  velocity  vg  is  also  aligned  with  the  grid 
elements  while  this  is  not  the  case  for  the  prismatic  grid.  Therefore,  if  the  preceding 
hypothesis  is  correct,  then  the  predictions  with  the  prismatic  grid  should  degrade  as  (3 
is  increased,  thus  increasing  the  relative  importance  of  vg.  That  this  is  the  case  can  be 
seen  from  Figs.  A.4-A.9  which  show  the  predictions  with  the  circulation  parameter 
(3  set  to  0.446  and  4.461.  As  (3  is  increased,  the  predictions  using  the  hexahedral 
grids  are  the  best,  then  those  of  the  hybrid  and  prismatic  grids,  then  those  of  the 
tetrahedral  grids.  The  tetrahedral  grids  degrade  the  least  as  (3  is  increased,  but  only 
because  their  predictions  were  so  poor  for  the  smallest  value  of  (3.  Fig.  8  also  shows 
the  relative  accuracy  of  the  four  grid  types  when  the  circulation  is  large.  It  compares 
the  pressure  and  circumferential  velocity  predicted  using  the  four  grid  types  with 
Nc  =  10  and  (3  =  4.461.  For  comparison,  the  prediction  using  the  hexahedral  grid 
with  Nc  =  22  is  also  shown;  it  is  close  to  the  true  solution. 

To  check  the  adequacy  of  the  resolution  in  the  x  direction  of  the  hexahedral,  hybrid 
and  prismatic  grids,  calculations  were  made  using  Nc  =  22  but  with  Nx  =  201  and 
401.  The  difference  in  the  calculated  pressures  and  velocities  is  barely  discernible 
when  plotted.  They  are  not  shown  here. 

The  calculations  on  the  tetrahedral  grid  are  consistent  with  an  artificial  viscosity 
caused  by  the  discretization  of  ^-momentum  of 

Ut 

^~0.01-A=  0.005C/A.  (11) 

where  A  =  2rc/Nc  is  the  size  of  the  cells  in  the  core.  For  an  accurate  solution 
the  artificial  viscosity  must  be  significantly  smaller  than  the  kinematic  viscosity,  u. 
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Therefore, 


O.OlUr, 


(12) 


V 


V 


A  <  200- 


which,  for  the  cases  described  here,  requires  that  Nc  22. 


When  (3  is  large,  the  discretization  of  circumferential  momentum  also  induces  sig¬ 
nificant  artificial  viscosity  on  the  tetrahedral,  prismatic  and  hybrid  grids.  In  that 
case: 


VQ, art 


0.01 


T c^O  max 

Nc 


0.001C 

Nc 


0.0005CA 

rc 


0.0005/3 f/A. 


(13) 


For  accurate  solutions  one  then  needs 


Nc  >  0.001 


c 

v 


0.001 


(3Urc 


A  < 


2000 vrc 

C 


2000z/ 

(3U 


(14) 


which  for  the  cases  described  here  requires  that  Nc  /$>  2.2/3. 


From  Eqs.  (11)  and  (13),  vxart  =  VQari  when  (3  —  10  which  suggests  that  the  hybrid, 
prismatic  and  tetrahedral  grids  should  have  similar  accuracy  when  f3  >  10,  but  that 
the  hybrid  and  prismatic  grids  will  be  more  accurate  for  smaller  (3. 


2.8  The  effects  of  the  discretization  scheme  for  the 
convective  terms 

The  calculations  shown  in  Annexes  A  and  B  were  calculated  using  the  ANSYS  CFX 
“High  Resolution  Scheme”  for  discretization  of  the  convection  terms  in  the  equations 
of  motion.  This  scheme  is  a  blend  of  a  second  order  accurate  scheme  and  a  first 
order  accurate  upwind  differencing  scheme.  The  blend  factor  is  calculated  locally  so 
the  scheme  is  as  close  to  second  order  accurate  as  possible  without  introducing  local 
oscillations  [6] . 

It  is  natural  to  assume  that  some  of  the  high  diffusivity  exhibited  by  the  solutions 
on  the  tetrahedral  grid  is  attributable  to  the  first  order  component  of  the  the  High 
Resolution  Scheme.  Therefore  all  the  laminar  calculations  on  the  tetrahedral  grid 
were  repeated  using  a  fully  second  order  accurate  scheme.  Some  of  the  calculations 
on  the  other  grids  were  also  repeated.  Fig.  9  compares  the  results  of  the  two  different 
discretization  schemes  for  the  case  with  (3  =  0.446  and  Nc  =  10  on  the  prismatic 
and  tetrahedral  grids.  The  results  on  the  prismatic  grid  show  very  little  change 
with  the  discretization  scheme:  so  much  so  that  it  is  often  difficult  to  discern  the 
symbols  for  the  High  Resolution  Scheme  as  they  are  overlaid  by  the  symbols  from 
the  second  order  accurate  scheme.  There  is  more  variation  with  the  tetrahedral 
grid  but,  counterintuitively,  the  second  order  scheme  increases  the  diffusivity  causing 
poorer  accuracy  relative  to  the  theoretical  curves.  Only  the  radial  component  of 
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Figure  9:  Comparison  of  the  High  Resolution  Scheme  and  a  fully  second  order  accu¬ 
rate  scheme  on  the  prismatic  and  tetrahedral  grids:  /3  =  0.446,  Nc  =  10.  The  top 
four  hgures  plot  pressure  coefficient,  circumferential  velocity,  axial  velocity  deficit  and 
radial  velocity  as  a  function  of  r/rCi  in  the  midstream  plane,  x  =  xmid .  The  bottom 
two  hgures  plot  pressure  coefficient  and  axial  velocity  dehcit  along  the  centreline. 
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the  velocity  shows  a  slight  improvement  with  less  scatter  when  the  High  resolution 
Scheme  is  used.  These  results  are  consistent  across  all  the  cases  tried:  the  hexahedral, 
prismatic  and  hybrid  grids  are  insensitive  to  the  discretization;  the  tetrahedral  grid 
is  more  sensitive  to  the  scheme  used  and  the  second  order  scheme  is  more  diffusive. 
We  do  not  understand  why  the  second  order  scheme  is  more  diffusive. 


3  The  evolution  of  a  turbulent  vortex 

3.1  RANS  calculations 

ANSYS  CFX  was  used  to  calculate  the  evolution  of  turbulent  vortices  in  water  a 
standard  temperature  and  pressure.  The  same  grids  were  used  as  in  the  laminar 
calculations.  The  upstream  plane  was  an  inlet  with  velocity  set  using  Eq.  (1)  with 
U  =  1  m/sec.  Three  values  for  the  circulation  were  used,  set  so  that  their  non- 
dimensional  values  were  the  same  as  for  the  laminar  case:  (3  =  0.0446,  0.446  and 
4.461.  The  outer  boundary  was  a  free  slip  wall.  The  downstream  plane  was  an  outlet 
with  average  normal  velocity  set  to  U.  In  the  upstream  plane  the  turbulent  kinetic 
energy  was  set  to 

k  =  |/f/2;  I  =  0.001  +  0.01e-(r/rci)2)  (15) 

and  the  eddy  viscosity  ratio  to 


^  =  1  +  1000e"(r/rci)2  (16) 

Since  the  default  kinematic  viscosity  of  water  used  by  ANSYS  CFX  is  8.926  x 
10-'  m2/sec,  the  maximum  value  of  Vt/UrCi  in  the  initial  plane  is  3.98  x  10”4,  close 
to  the  value  of  v/UrCi  for  the  laminar  runs  (4.46  x  10~4).  However,  in  the  turbulent 
case  ut  is  confined  within  the  vortex  core  and  evolves  with  the  flow. 

The  Shear  Stress  Transport  (SST)  turbulence  model  was  used  for  all  turbulent  cal¬ 
culations. 

3.2  Results  of  calculations 

Annex  B  contains  figures  showing  the  results  of  the  turbulent  calculations  on  all  the 
grids  for  flows  with  (3  =  0.0446,  0.446  and  4.461.  The  predicted  values  of  pressure 
coefficient,  Cp,  circumferential  velocity,  vq/U,  axial  velocity  deficit,  (U  —  vx)/U,  and 
vt/Urci)  are  plotted  as  a  function  of  r/ra  in  the  plane  x  =  xmui.  The  pressure  and 
turbulent  viscosity  are  also  plotted  as  a  function  of  (x  —  x0)/rci  along  the  centreline 
(r  =  0). 
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Figure  10:  The  variation  of  the  turbulent  kinematic  viscosity  along  the  cen¬ 
treline  for  different  values  of  the  circulation.  On  the  left  the  data  is  plotted 
as  vt  vs.  (x  —  xq )/rCi;  on  the  right  it  is  plotted  as  vt/C  vs.  (3{x  —  xq )/rCi. 

The  only  significant  difference  in  the  way  the  laminar  and  viscous  calculations  were 
set  up  is  in  the  kinematic  viscosity.  Therefore  we  first  consider  how  it  changes  as  the 
vortex  evolves. 

In  a  turbulent  vortex,  if  the  turbulent  viscosity  is  significantly  higher  than  the  laminar 
viscosity,  then  in  the  cross-stream  plane,  in  which  the  flow  is  largely  decoupled  from 
the  flow  in  the  downstream  direction,  there  are  only  two  dimensional  parameters 
of  significance:  the  circulation,  C,  and  the  core  radius,  rc.  This  implies  that  the 
turbulent  viscosity  will  be  roughly  proportional  to  C:  vt~  aC .  If  vt  it  is  too  high,  the 
turbulent  viscosity  diffuses  radially,  smearing  itself  over  a  larger  area  until  its  diffusion 
is  just  matched  by  its  production;  if  it  is  too  low,  production  will  dominate  diffusion, 
increasing  the  turbulent  viscosity  until  diffusion  and  production  are  balanced. 

In  practice,  when  the  turbulent  viscosity  is  initially  smaller  than  aC,  it  may  take 
a  long  time  before  it  increases  to  the  point  where  diffusion  balances  production. 
The  time  scale  over  which  it  increases  is  proportional  to  the  inverse  of  the  mean 
strain.  Since  the  mean  strain  is  proportional  to  C/r^  =  /3U/rc,  the  time  scale  is 
inversely  proportional  to  /3;  weaker  vortices  take  longer  to  evolve.  Put  another  way, 
the  distance  traversed  before  the  core  begins  its  rapid  production  of  turbulent  viscosity 
is  proportional  to  (3. 

The  left  side  of  Fig.  10  shows  the  turbulent  kinematic  viscosity  along  the  centreline  for 
three  different  values  of  (3  on  the  hexahedral  grid  with  22  nodes  across  the  core  (extra 
calculations  were  done  to  include  the  intermediate  (3  values  of  0.1338  and  1.338).  The 
calculations  were  performed  using  the  hexahedral  grid  with  22  nodes  across  the  core. 
When  j3  =  4.461,  after  a  small  initial  decrease,  the  turbulent  viscosity  increases  very 
rapidly  at  about  50  core  radii  downstream  of  the  initial  plane;  for  [3  =  1.338  the  rapid 
increase  is  delayed  until  about  150  core  radii;  for  f3  =  0.446  the  viscosity  is  only  just 
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beginning  to  increase  at  450  core  radii.  At  the  two  lower  values  of  [3  the  viscosity 
continues  to  decrease  for  450  core  radii. 

The  right  side  of  Fig.  10  shows  the  same  data  as  the  left  side  but  with  vt/C  plotted 
against  (3{x  —  Xo)/rCi.  (Additional  calculations  were  required  to  make  this  figure, 
both  at  the  intermediate  values  of  (3  and  on  a  grid  elongated  to  x\  =  5000  metres. 
Even  so  the  curves  at  the  lower  values  of  f3  are  truncated  because  the  grid  was  not 
long  enough.)  Despite  the  differing  values  of  vt/C  at  the  initial  plane,  for  each  value 
of  C  the  rapid  production  of  vorticity  begins  near  (3(x  —  x0 )/rci  =  200.  Moreover, 
after  the  rapid  production  the  value  of  vt j C  appear  to  be  approaching  a  constant 
value,  a,  near  0.004,  though  this  is  uncertain  for  the  lower  values  of  (3  for  which  the 
curves  are  truncated. 


The  value  of  a  at  which  production  matches  diffusion  is,  in  the  computer  codes, 
dependent  on  the  turbulence  model  used,  ft  is  well  known  that  isotropic  turbulence 
models  tend  to  overpredict  the  production  of  turbulent  viscosity  in  vortices;  therefore 
the  value  of  0.004  suggested  by  these  results  for  the  SST  model  is  expected  to  be 
higher  than  the  values  predicted  by  non-isotropic  models. 

Once  diffusion  and  production  are  matched,  the  viscosity  is  roughly  constant  and  we 
can  expect  it  to  evolve  roughly  as  a  laminar  vortex.  Therefore,  from  Eq.  (7),  its  core 
radius  will  double  in  a  distance 

r  1.6  Ur?  1.6  rc 

vt  <*(3 

and  over  the  same  distance  the  maximum  velocity  will  reduce  by  half  and  the  pressure 
drop  will  decrease  to  one  quarter  of  its  value.  For  values  of  [3  near  1.0,  Lr  is  a  few 
hundred  core  radii.  For  smaller  values  of  (3  it  can  be  thousands  of  core  radii. 


For  the  lowest  value  of  / 3  in  the  numerical  calculations,  the  initial  value  of  vt  exceeds 
aC.  Therefore  the  viscosity  is  reduced  by  self-diffusion  resulting  in  a  vortex  whose 
effective  Reynolds  number,  Urjvf,  is  larger  than  for  the  laminar  case.  Therefore  we 
can  expect  the  results  on  the  tetrahedral  grid  to  be  very  poor,  the  turbulent  viscosity 
being  completely  masked  by  artificial  viscosity.  Figs.  B.1-B.3  show  this  to  be  true. 
Moreover,  as  can  be  seen  in  the  lower  right  graphs  of  Figs.  B.2  and  B.3,  the  artificial 
viscosity  diffuses  the  turbulent  viscosity  so  much  that  it  is  effectively  removed.  On  the 
other  hand,  the  other  three  grids  all  give  very  accurate  results  for  this  case,  converging 
well,  as  in  the  laminar  case,  with  only  a  few  nodes  across  the  core. 

When  the  strength  of  the  vortex  is  increased  to  (3  =  0.446  (Figs.  B.4-B.6),  the  initial 
value  of  vt  is  significantly  smaller  aC,  so  that  it  will  eventually  increase.  From  Fig.  B.6 
we  can  see  that  at  the  mid-plane  the  increase  is  just  beginning.  The  graphs  on  the 
right  of  Fig.  B.5  show  that  viscosity  is  being  produced  near  the  location  of  peak 
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velocity  where  the  mean  strain  is  highest.  For  this  value  of  (3  the  trends  with  varying 
number  of  nodes  are  again  similar  to  the  laminar  case:  the  hexahedral  grid  converges 
very  well,  the  hybrid  and  prismatic  grids  are  similar  but  not  quite  as  good  as  the 
hexahedral  grid,  and  the  tetrahedral  grid  is  still  very  poor. 

When  the  strength  of  the  vortex  is  increased  further  to  3  —  4.461  (Figs.  B.7-B.9)  the 
flow  in  the  mid-plane  is  significantly  different  from  the  cases  with  lower  (3 .  Now  the 
vorticity  has  already  undergone  the  major  portion  of  its  increase  and  is  approaching 
its  final  state  as  a  well-developed  vortex.  An  interesting  feature  of  the  flow  is  that 
increasing  the  node  density  now  results  in  lower  maximum  velocities  and  smaller 
pressure  drops  in  the  core.  The  reason  for  this  can  be  found  in  Fig.  B.9.  When  the 
node  density  is  lower,  increased  artificial  viscosity  results  in  greater  diffusion  of  the 
velocity  in  the  early  stages  of  evolution.  The  mean  strain  is  then  lower  and  the  rapid 
increase  of  vorticity  is  delayed.  As  a  result  the  vortex  does  not  diffuse  as  much  and 
the  maximum  velocity  is  higher. 

Judging  the  convergence  characteristics  from  Figs.  B.7-B.8,  it  is  difficult  to  choose 
between  the  three  non-tetrahedral  grids,  but  Fig.  B.9  once  again  shows  that  the 
hexahedral  grid  is  best  and  the  hybrid  and  prismatic  grids  are  somewhat  poorer.  As 
always  the  tetrahedral  grid  ranks  last. 

4  Implications  for  the  gridding  of  propellers 


The  flow  past  propeller  DTMB  4119,  shown  in  Fig.  11,  has  been  studied  experimen¬ 
tally  by  Jessup  et  al.  [7,8]  and  was  used  as  a  test  case  by  the  22nd  ITTC  Propulsion 
Committee  Propeller  RANS/Panel  Method  Workshop  in  Grenoble,  1998. 

Let  V  be  the  speed  of  advance  of  the  propeller  and  let  Vr  be  the  component  of  the 
fluid  velocity  pointing  radially  outward  from  the  propeller  axis.  Fig.  12  plots  measured 
values  of  Vr/V  for  DTMB  4119  on  a  circle  just  downstream  of  the  propeller  blades, 
centred  on  the  propeller  axis  and  passing  through  the  core  of  the  tip  vortices.  The 
location  of  the  circle  is  shown  in  Fig.  11.  As  this  is  a  three-bladed  propeller,  Vr  is 
only  plotted  along  one  third  of  the  circle.  The  core  of  the  vortex  is  easily  identified  as 
the  region  between  the  maximum  and  minimum  values  of  Fn  a  region  with  angular 
extent  of  about  4  degrees.  As  the  circle  has  a  radius  of  r  =  0.9247?,  this  represents  a 
distance  of  about  0.0657?,  where  7?  is  the  propeller  radius.  However,  the  vortex  cuts 
through  the  plane  of  the  circle  obliquely  causing  the  vortex  to  be  stretched  by  the 
factor  \Jl  +  (nr/ JR)2  in  the  direction  tangent  to  the  circle.  Since  the  advance  ratio 
J  =  0.833,  this  means  that  the  actual  diameter  of  the  core  is  about  3.6  times  smaller 
than  the  distance  on  the  plot:  i.e.  rc  ~  0.0097?. 

The  vortex  strength,  /?,  can  also  be  estimated  from  the  difference  in  value  of  the 
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Figure  11:  The  propeller  DTMB  4119:  looking  directly  upstream 
(left)  and  from  the  side  (right).  The  yellow  circle  marks  the  location 
of  the  data  in  Figure  12. 


Figure  12:  The  radial  component  (relative  to  the  propeller  axis)  of 
the  flow  past  DTMB  4119  along  a  circle  passing  through  the  vortex 
core. 
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minimum  and  maximum  of  Vr  and  Eq.  (6).  Note  that  the  convection  speed,  U,  is 
the  combination  of  the  inflow  speed,  V,  and  the  velocity  induced  by  the  propeller 
rotation:  U  =  Vyjl  +  (nr/ JR)2  «  3.6V.  Therefore 


A 14 


Vdi 


0.776E;  (5  = 


C 

U^r. 


Voi 


0.1147/ 


V 

6.8—  «  1.9 


(18) 


Since  V  =  2.54  m/sec  and  the  propeller  radius  is  152.4  mm,  the  circulation  can 
be  estimated  as  C  —  2.4  x  10-2  m2/sec  and  the  limiting  value  of  the  turbulent 
viscosity  once  the  vortex  is  well-developed  will  be  about  ut  ~  0.004(7  m  10-4  m2/sec. 
Thereafter  the  maximum  radius  will  decrease  by  half  over  a  distance  of  roughly  2 R 
or  about  one  third  of  a  turn  of  the  propeller. 


Fig.  10  suggests  that  the  vortex  will  evolve  slowly  for  about  100  core  radii  before 
the  vorticity  increases  sufficiently  that  it  diffuses  quickly.  What  is  the  best  way  to 
construct  the  grid  if  we  wish  to  predict  the  vortex  characteristics  accurately  over  this 
distance? 


A  tetrahedral  grid  is  easiest  to  generate,  but  it  will  require  a  minimum  of  20  nodes 
across  the  core  which  gives  about  600,000  nodes  just  within  the  core  of  a  single  vortex. 
Moreover,  DTMB  4119  is  a  model  scale  propeller;  for  propellers  at  full  scale  the  vortex 
core  is  expected  to  be  even  smaller.  Therefore,  even  if  grid  adaptation  is  used  to  refine 
the  grid  near  the  vortex  cores,  very  large  grids  will  be  required. 

Although  the  numerical  tests  suggest  that  a  hexahedral  grid  yields  the  best  accuracy 
for  the  fewest  number  of  nodes,  that  is  really  only  because  the  grid  is  perfectly  centred 
on  the  vortex.  As  soon  as  the  vortex  core  moves  away  from  the  centre  of  the  grid, 
the  velocity  components  cross  the  element  faces  obliquely  and  the  hexahedral  grids 
give  results  similar  to  the  prismatic  grid.  Therefore  the  location  of  the  vortex  would 
meeds  to  be  known  very  accurately  if  the  benefits  of  the  hexahedral  grid  are  to  be 
obtained.  Moreover,  the  upstream  and  downstream  faces  of  the  hexahedral  grids 
contain  very  thin  triangular  elements.  These  make  it  hard  to  embed  a  hexahedral 
grid  within  a  larger  unstructured  grid  covering  the  remainder  of  the  fluid.  Therefore, 
the  hexahedral  grid  is  not  a  practical  option. 

On  the  other  hand,  only  a  reasonable  estimate  of  the  locations  of  the  vortex  cores  is 
necessary  to  get  adequate  alignment  of  the  hybrid  or  prismatic  grids.  Such  estimates 
are  easily  obtained  from  actuator  disk  theory  once  the  thrust  is  known.  These  grids 
offer  the  prospect  of  a  considerable  reduction  in  the  number  of  nodes,  both  because 
somewhat  fewer  nodes  are  required  across  the  core  to  get  similar  accuracy,  and  because 
the  elements  can  be  elongated  along  the  direction  of  the  vortex.  Moreover,  there  is 
no  difficulty  in  embedding  a  hybrid  or  prismatic  grid  enclosing  the  vortex  within  a 
larger  unstructured  grid.  Therefore  these  grid  types  are  the  most  attractive  option.  A 
simple  gridding  procedure  would  be  to  calculate  the  flow  on  an  initial  grid  to  estimate 
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the  thrust,  then  generate  a  new  grid  with  prismatic  refinement  along  the  vortex  cores. 
One  can  easily  envisage  further  iterations  of  the  grid  as  the  vortex  locations  become 
known  more  accurately  with  each  set  of  calculations. 


5  Concluding  remarks 


The  evolution  of  laminar  and  turbulent  rectilinear  vortices  has  been  studied  using  the 
flow  solver  ANSYS  CFX.  The  results  confirm  the  finding  of  an  earlier  study  by  Hally 
and  Watt  [1]  that  accurate  predictions  of  weak  laminar  vortices  (/?  -C  1)  over  distances 
equal  to  hundreds  of  core  radii  can  be  attained  with  6  to  10  nodes  across  the  core. 
However,  good  predictions  can  only  be  obtained  if  the  grid  is  comprised  of  hexahedral 
or  prismatic  elements  whose  faces  are  either  aligned  with,  or  perpendicular  to,  the 
vortex  axis.  In  cases  in  which  the  faces  cut  the  vortex  at  large  oblique  angles,  in 
particular  a  tetrahedral  grid,  artificial  dissipation  induced  by  the  convection  velocity 
across  the  faces  results  in  excessive  dissipation  and  poor  predictions. 

The  current  study  also  examined  strong  laminar  vortices  with  /3  of  the  order  of  1 
or  greater.  For  many  practical  flows,  including  propeller  tip  vortices,  this  is  a  more 
realistic  range  of  values.  In  a  strong  vortex  the  circumferential  velocity  is  similar 
in  magnitude  to  the  convection  velocity  and  the  artificial  dissipation  caused  by  the 
circumferential  velocity  crossing  oblique  faces  in  the  grid  becomes  an  important  fac¬ 
tor.  Therefore,  even  when  the  element  faces  are  aligned  with  the  vortex  axis,  larger 
numbers  of  nodes  are  needed  to  reduce  the  artificial  dissipation  to  acceptable  levels: 
when  f3  =  0.446,  approximately  16  nodes  across  the  core  and,  when  f3  =  4.446,  more 
than  22  nodes  across  the  core  to  get  similar  accuracy  as  6  nodes  across  the  core  when 
(3  =  0.0446.  Nevertheless,  one  still  finds  that  there  is  a  significant  improvement  when 
the  grids  are  aligned  (the  hexahedral,  hybrid  and  prismatic  grids)  than  when  they 
are  not  (the  tetrahedral  grid). 

When  the  vortex  is  turbulent,  it  will  evolve  towards  a  well-developed  state  in  which 
the  turbulent  vorticity  is  proportional  to  the  circulation.  Therefore  weak  vortices  will 
have  low  levels  of  vorticity  and  will  consequently  evolve  slowly.  The  effects  of  the  grid 
on  the  flow  predictions  are  similar  to  those  of  the  laminar  case.  When  the  vortex  is 
strong,  an  initial  low  level  of  vorticity  will  increase  towards  its  final  state;  the  stronger 
the  vortex,  the  more  rapidly  the  vortex  will  reach  its  fully-developed  state.  In  these 
cases  artificial  viscosity  can  have  an  important  effect  in  reducing  the  mean  strain  that 
drives  the  production  of  vorticity;  therefore,  coarse  grids  overestimate  the  time  taken 
to  achieve  the  well-developed  state. 
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Annex  A:  Results  of  the  laminar  computations 


The  figures  on  the  following  pages  show  the  results  of  the  laminar  computations  for 
each  type  of  grid.  For  each  grid  and  for  each  value  of  the  parameter  /?,  six  plots  are 
presented: 

1.  the  pressure  coefficient  as  a  function  of  r/rci  in  the  plane  x  =  xmid, 

2.  the  circumferential  velocity,  Ve/U,  as  a  function  of  r/rci  in  the  plane  x  =  xnuci] 

3.  the  axial  velocity  deficit,  (U—vx)/U,  as  a  function  of  r/ra  in  the  plane  x  =  xmid ; 

4.  the  radial  velocity,  vr/U,  as  a  function  of  r/rCi  in  the  plane  x  =  xmid ; 

5.  the  pressure  coefficient  along  the  centreline  as  a  function  of  the  distance  from 
the  initial  plane  measured  in  core  radii  (i.e.  (x  —  x0 )/rci);  and 

6.  the  axial  velocity  deficit,  (U  —  vx)/U,  along  the  centreline  as  a  function  of  the 
distance  from  the  initial  plane  measured  in  core  radii. 

Each  plot  shows  the  results  for  six  runs  with  varying  numbers  of  nodes  across  the 
vortex  core. 

For  easy  comparison  between  the  different  grids,  the  results  for  each  of  the  three  grid 
types  are  shown  on  the  same  page. 

In  all  the  figures,  the  solid  black  line  shows  the  values  predicted  using  Eq.  (1).  When 
/3  =  0.0446  this  provides  an  accurate  estimate  of  the  true  solution;  all  numerical 
calculations  should  converge  to  the  black  line.  When  f3  =  0.446,  the  approximations 
used  in  Eq.  (1)  are  at  the  limits  of  their  validity;  all  numerical  calculations  should 
converge  close  to  the  black  line  but  small  discrepancies  should  be  expected  within 
the  vortex  core.  When  (3  =  4.461,  the  approximations  used  in  Eq.  (1)  are  no  longer 
valid;  one  should  not  expect  the  numerical  predictions  to  converge  to  the  black  line; 
however,  it  is  retained  in  the  plots  to  make  it  easier  to  compare  between  the  results 
at  different  circulations. 
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Figure  A.2:  Laminar  vortex:  the  axial  velocity  deficit  (left)  and  radial  velocity 
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Figure  A.8:  Laminar  vortex:  the  axial  velocity  deficit  (left)  and  radial  velocity 
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Annex  B:  Results  of  the  turbulent 
computations 


The  figures  on  the  following  pages  show  the  results  of  the  turbulent  computations  for 
each  type  of  grid.  For  each  grid  and  for  each  value  of  the  circulation  parameter  /3, 
six  plots  are  presented: 

1.  pressure  coefficient  as  a  function  of  r/rci  in  the  plane  x  =  xmid ; 

2.  the  circumferential  velocity,  vq/U,  as  a  function  of  r /rcl  in  the  plane  x  =  xmid ; 

3.  the  axial  velocity  deficit,  (U—vx)/U,  as  a  function  of  r/rCj  in  the  plane  x  =  xmuj; 

4.  the  non-dimensional  turbulent  viscosity,  vt / Urci ,  as  a  function  of  r/rCj  in  the 
plane  x  =  xmid  metres; 

5.  the  pressure  coefficient  along  the  centreline  as  a  function  of  the  distance  from 
the  initial  plane  measured  in  core  radii  (i.e.  (x  —  x0 )/rci);  and 

6.  the  non-dimensional  turbulent  viscosity,  vt/Ur ci,  along  the  centreline,  as  a  func¬ 
tion  of  the  distance  from  the  initial  plane  measured  in  core  radii. 

Each  plot  shows  the  results  for  six  runs  with  varying  numbers  of  nodes  across  the 
vortex  core. 

For  easy  comparison  between  the  different  grids,  the  results  for  each  of  the  three  grid 
types  are  shown  on  the  same  page. 
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Figure  B.1 :  Turbulent  vortex:  the  pressure  coefficient  (left)  and  circumferential 
velocity  (right)  as  a  function  of  r/rCi  in  the  plane  x  =  xmid  with  (3  =  0.0446. 
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Figure  B.2:  Turbulent  vortex:  the  axial  velocity  deficit  (left)  and  the  turbulent 
viscosity  (right)  as  a  function  of  r/rCi  in  the  plane  x  =  xmid  with  j3  =  0.0446. 
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Figure  B.3:  Turbulent  vortex:  the  pressure  coefficient  (left)  and  the  turbulent 
viscosity  (right.)  as  a  function  of  (x  —  Xo)/rCi  on  the  centreline:  f3  =  0.0446. 
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Figure  B. 4:  Turbulent  vortex:  the  pressure  coefficient  (left)  and  circumferential 
velocity  vq/U  (right)  as  a  function  of  r/rci  in  the  plane  x  =  xmid  with  (3  =  0.446. 


DRDC  Atlantic  TM  2009-052 


37 


n/(xA-n)  n/(xA-n)  n/(xA-n)  n/(xA-n) 


Hexahedral  Grid:  p  =  0.446 


Hexahedral  Grid:  p  =  0.446 


0.0014 

0.0012 

0.001 

0.0008 

0.0006 

0.0004 

0.0002 

0 

-0.0002 


0.0003 

3  nodes 

A 

3  nodes 

6  nodes 

0.00025 

6  nodes 

1 0  nodes  ** 

j 

1 0  nodes  « 

1 6  nodes  ° 

16  nodes  » 

22  nodes 

0.0002 

6 

22  nodes 

L  x 

0.00015 

\ 

> 

o 

0.0001 

-  4  + 

6 

5e-05 

A 

, 

0 

% 

. 

0123456789 


0123456789 


r/rci 


Hybrid  Grid:  p  =  0.446 


Hybrid  Grid:  p  =  0.446 


4  5 

r/iri 


Prismatic  Grid:  p  =  0.446 


0123456789 


r/rd 


Tetrahedral  Grid:  p  =  0.446 


Tetrahedral  Grid:  p  =  0.446 


0.0003 

- 1 - 1 - 1 - , - [ - [ - [ - , - 

3  nodes 

0.00025 

6  nodes 

1 0  nodes  » 

16  nodes  ° 

0.0002 

22  nodes 

0.00015 

0.0001 

5e-05 

0 

)1  23456789 

Figure  B.5:  Turbulent  vortex:  the  axial  velocity  dehcit  (left)  and  the  turbulent 
viscosity  (right)  as  a  function  of  r/rCi  in  the  plane  x  =  xmid  with  j3  =  0.446. 
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Figure  B.6:  Turbulent  vortex:  the  pressure  coefficient  (left)  and  the  turbulent 
viscosity  (right)  as  a  function  of  ( x  —  Xq )/rCi  on  the  centreline:  [3  =  0.446. 
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Figure  B.  7:  Turbulent  vortex:  the  pressure  coefficient  (left)  and  circumferential 
velocity  (right)  as  a  function  of  r/rCi  in  the  plane  x  =  xmid  with  /3  =  4.461. 
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Figure  B.8:  Turbulent  vortex:  the  axial  velocity  de&cit  (left)  and  the  turbulent 
viscosity  (right)  as  a  function  of  r/rCi  in  the  plane  x  =  xmid  with  j3  =  4.461. 
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Figure  B.9:  Turbulent  vortex:  the  pressure  coefficient  (left)  and  the  turbulent 
viscosity  (right)  as  a  function  of  ( x  —  Xq )/rCi  on  the  centreline:  [3  =  4.461. 
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Annex  C:  Perturbation  expansion  of  the 
vortex  solution 


To  get  a  firm  grasp  on  the  size  of  the  errors  in  the  approximate  vortex  solution,  we 
express  this  solution  as  a  perturbation  expansion  in  a  small  parameter  e.  We  assume 
that  U  is  of  order  1  and  vg  is  of  order  e.  From  Eq.  (6),  that  implies  that 

V°^-  =  0.114—  =  0.114/3  (C.l) 

U  Urc 

is  of  order  e.  Therefore  the  approximation  is  valid  if  (3  <C  9.  If  we  let  rc  be  of  order 
1  (this  simply  sets  a  length  scale  for  the  flow),  then  C  must  be  of  order  e. 

Furthermore,  to  model  the  dominance  of  convection  over  diffusion,  we  assume  that 
derivatives  with  respect  to  x  are  of  order  e:  i.e.  rc/x  is  of  order  e.  Since  rc  ~ 
2.24a/ vx/U,  this  implies  that  5u/Urc  -C  1  and  v  is  of  order  e. 

Let 


C  =  eC*;  v  =  eu*- 

x  =  x* /e 

(C.2) 

and  assume  a  solution  of  the  form: 

vx{x*,r)  = 

\x*,r) 

(C.3) 

vr(x*,r)  =  Y  ezv^ 

\x*,r) 

(C.4) 

vg(x*,r )  = 

\x*,r) 

(C.5) 

p(x*,r )  =  ^  elp^1 

i 

\x*,r) 

(C.6) 

To  conform  with  the  approximate  solution  we  have 

v{x\x*,r)  =  U 

(C.7) 

v^\x*,r)  =  v^\x*,r)  =  0 

(C.8) 

O 

* 

o 

(C.9) 

„<■> = f  (l  - 
6  2vr r  V 

Ur2  \ 
e  x*  ) 

(C.10) 

(x* ,  r)  =  p (x*,r)  =  0 

(C.ll) 

p,2>(x'’r)-  4?r2  l 

00  1  /  u *2  \2 

— ^  /I  —  e  i"***  J  dz 

(C.12) 

From  here  on  we  will  drop  the  asterisks  on  the  perturbed  variables.  Since  we  will  no 
longer  use  the  original  variables,  this  should  cause  no  confusion. 
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The  Navier-Stokes  equations  for  steady  axisymmetric  flow  in  cylindrical  coordinates 
can  be  written: 


dvx  1  d(rvr ) 
dx  r  dr 

dvx  dvx 
—  +  ry— 
ox  or 


dvr 


dvr 


dx  Vr  dr 

dvg  dvg 

evx  —  +  vr—  + 
dx  dr 


_  VA 

r 

vrve 

r 


0 


_e  dp  f  2d2vx 

p  dx  6  \  dx2 

1  dp  (  2  d2vr 

pdr  eV  \e  dx2 

/A  in 

\  dx2  r  dr  \ 


1  d_ 

r  dr 


1  d_ 

r  dr 


(C.13) 

(C.14) 

(C.15) 

(C.16) 


Substituting  Eqs.  (C.7)-(C.12)  into  Eqs.  (C.13)-(C.16)  and  retaining  only  the  lowest 
order  of  e  in  each  equation  we  get: 


dv^ 

dx 


1  d(rvr 
r  dr 


(2)> 


u 


dv. 


(i) 


dx 


0 

v  d  (  dv^ 

- r - 

r  dr  \  dr 


(i^)2  1  dp^ 

r  p  dr 

Ud-f  =  v(~ 

dx  \  r  dr 


dv. 


(i) 


dr 


(C.17) 

(C.18) 

(C.19) 

(C.20) 


Eqs.  (C.19)  and  (C.20)  have  solutions  given  by  Eqs.  (C.10)  and  (C.12).  Eq.  (C.18)  can 
be  solved  for  vx  ^ ;  the  solution  with  vx1\xo,r)  =  0  is  simply  vx1\x,r)  =  0.  Eq.  (C.17) 
then  implies  that  vx2\x,r )  =  0. 


Proceeding  to  the  next  level  of  expansion  we  get: 

dv i2'*  1  d(ri’r^) 

dx  r  dr 

TTdv 1  dpt2'1  v  d  (  dv i2'* 

U -  — - — (— - I  T“ - 

dx  p  dx  r  dr  y  dr 

2v^Vg2)  1  dp ^ 

r  p  dr 

dv 02)  ( 1  d  (  dvf} 

U~d^T  =  U  [  ~r~dr  |r^T 


(C.21) 

(C.22) 

(C.23) 

(C.24) 
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/q\  /"oN 

Eq.  (C.24)  implies  that  vye  =  0  to  satisfy  the  initial  condition  that  v#  (xo ,  r)  =  0. 
Then  Eq.  (C.23)  implies  that  =  0. 


Let 


Then 


where 


= 


Pc2  r°°  i 

47T2  Jr  z3 
/(*)  = 


r2U 
z  =  A 

Aux 

(C.25) 

^1  —  e~  4^^  dz  =  - 

^2f/  ft  ^ 
‘32IrVx/(") 

(C.26) 

riy?)'* 

(C.27) 

(2) 

We  look  for  a  full  solution  to  Vx  ;  of  the  form 


32v t2u 

Substituting  into  Eq.  (C.22)  we  find  that 
(x  -  x0) 


^M  =  ^-±^Lan(z) 


71=1 


X 


(C.28) 


E 

n= 0 


X" 


((+  +  l)9n+i{z)  -  (n+l)gn(z)  -  (l  +  z)g'n(z)-zg"(z)'j 


provided  that  we  define  go{z)  =  0.  Therefore, 


9o(z)  =  0 

9i(z)  =  -/(-)  -  zf\z ) 

(1  +  z)g'n{z)  +  zg”(z) 


gn+ i{z)  =  gn(z 


71+1 


-m-zf(z) 

(C.29) 

(C.30) 

(C.31) 

(C.32) 


Let  hn(z)  =  ( gn{z )  -  gi(z))e2z.  Then 
C 2  ZZ  (x  —  x0)n 


327 t2u 

C 2 

327 rV 


n=  1 


X 


n+1 


z  e 


-2z 


) 


gi(t)(x-x„)  :  ,  _2_. 

+  ~n+l  an\^)e 


XXq 


n= 2 


C2(a;  -  go) 
327r2z/a;a:o 


/(*) 


(1  —  e  z)2  a:0e 


-22 


a: 


E 

n=  1 


X  —  Xq 


X 


hn+i  (z) 


(C.33) 
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The  functions  hn(z)  are  polynomials  in  z  which  can  be  determined  from  a  recursion 
obtained  from  the  recursion  relation  for  gn(z): 


hi(z)  =  0 

2(z  —  l)hn(z)  +  (1  —  3z)h'n(z)  +  zh'n(z)  +  2 

nn+l\Z)  ~  nn\Z)  -\ - — - 

71+1 

The  first  few  hn(z)  are: 


hi{z) 

h2(z) 

h3(z) 
h4(z) 
h5(z) 
K(z ) 
h7(z) 
h%{z) 


0 

A 

7  z  z^ 

6 +  3 +  y 

7  llz  _  z2  2z3 
6  +  iy  _  15  +  15" 

37  7z  7z2  _  2z3  2 z4 

3o  +  i5  +  iy_iy+4y 

37  79 z  llz2  14z3  26z4  4z5 

30  +  105  _  "l05~  +  ~45~  ~  "3I5"  +  315 
533  113z  113z2  97z3  113z4  llz5  z6 

420  +  ~210~  +  210  _  ~315  +  630  ~~  "315"  +  315 


(C.34) 

(C.35) 


(C.36) 

(C.37) 

(C.38) 

(C.39) 

(C.40) 

(C.41) 

(C.42) 

(C.43) 


Along  the  centreline  we  have 

C2 


.(2)/ 


vp(x,  0) 


32t r2 


v 


21 

\Xo  x  J  X4 

1  |  7(x  —  xp)2  |  37(x  —  xp)4  |  533(3;  -Xp)6  | 


6x2 


30x4 


420x6 


(C.44) 


/ON 

Eq.  (C.21)  can  now  be  used  to  solve  for  ry  . 

1  d(rvi^)  dvf^ 

r  dr  dx 

(o\ 

Use  Eq.  (C.22)  to  eliminate  the  term  in  dvyx  /dx, 


so  that 


1  dp^  v  d  f  dv ^ 
pU  dx  Ur  dr  l  dr 


1  f  dp  2)  v  dvi 

-  /  r  — ; —  dr - - — 

pUr  J o  dx  U  dr 


(C.45) 


(C.46) 


(C.47) 
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Now 


Therefore 


PC2U  d 

(M' 

32ti2u  dx 

V  x 

dpi2) 

— —  dr  = 

c 2 

327 t2isx2 


(C.48) 


r</< 

lQ7T2xrU  J o  ' 

c^m 

16n2xrU 
C2r  fr2U 
647 t2isx2  \  4 ux 


(f(z)  +  zf\z))rdr 
(f(z)  +  zf(z))  dz 


(C.49) 


Moreover 

dv^  dvx ?  dz  2z  dv^ 

dr  dz  dr  r  dz 

-  02 z  9  \(  f(~)  - 

167 T2vr  dz  I 


0--e~zf 

z 


1  1 

Xq  X 


n= 2 


C2z  2(1  —  e  z)e  z(a;  —  x0) 
167r2ixrx  zxn 


71=2 


where 

The  first  few  gn(z)  are 


9n(*)  =  ^n(-)  -  2hn(z) 


qi(z)  =  0 
q2(z)  =  -2 

4  4z 
ysv  7  33 

2^2 

94  (*)  =  -2  -  — 

,  v  8  8z  8z2  4z3 

^5V  ;  55  15  15 

.  .  0  4^2  4z3  4 z4 

^bV  ;  3  9  45 

,  ,  12  12z  8^2  20z3  8z4  8 z5 

,7W  =  _y_^r+^“_“2r+ 35  “sis 

,  ,  „  „  o  4^3  8z4  4^5  2z6 

(?8(,)  =  -2_2,2  +  ---  +  --- 


(C.50) 

(C.51) 

(C.52) 

(C.53) 

(C.54) 

(C.55) 

(C.56) 

(C.57) 

(C.58) 

(C.59) 
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Therefore, 


v3>  = 


C2r 


64ir2isx2 


/(*)" 


2(1  —  e  2)e  z(x  —  x0) 


zx  o 


E 

n= 2 


(x  -  x0y 


xri 


Qn(z)e 


—2z 


(C.60) 


(o\ 

Notice  that  tv  '  is  not  identically  zero  on  the  initial  plane.  This  is  because  the 
momentum  flux  due  to  the  gradient  of  p  in  the  axial  direction  must  be  balanced  by 
a  gradient  in  vx:  i.e.  dv^ /dx  ^  0  in  the  initial  plane.  The  continuity  equation  then 
requires  that  vr  ^  0  in  the  initial  plane. 


The  leading  terms  of  each  of  the  independent  variables  are  now  known. 

vx  =  JJ  +  v^e2  +  o(e3) 
vr  =  vl3)e3  +  o(e4) 
ve  =  evf]  +  o(e3) 
p  =  p^e2  +  o(e4) 


(C.61) 

(C.62) 

(C.63) 

(C.64) 


We  proceed  to  the  next  level  of  expansion  to  get  an  idea  of  the  error  in  these  terms. 

,(4h 


U- 


dvi 3)  ^  1  d/rv/1 
dx 


_| - 

r 


dr 


=  0 


dx 


1  <9p(4)  2v^vf')  TT<9vr3^ 


p  dr 


U- 


dx 


,  1  d  (  dv 
+  v  r 


(3)' 


(3)' 

Vr 


r  dr 


dr 


Lrdvf 

dx 


,(3) 


+  V 


=  V 


(2)dv(e]  (z)dv{e] 
dx2  x  dx  r  dr 


(C.65) 

(C.66) 

(C.67) 

(C.68) 


Eq.  (C.66)  implies  that  v ^  =  0,  since  it  is  zero  on  the  upstream  plane.  Eq.  (C.65) 
then  implies  that  v(r4'1  =  0.  Eq.  (C.68)  can  be  solved  using  the  known  expressions  for 
the  terms  on  the  right  hand  side.  Eq.  (C.67)  can  then  also  be  solved. 


The  first  term  on  the  right  hand  side  of  Eq.  (C.68)  is 

Ur 2 


d2v ^ 

/ - - - 

dx2 


CUr 


327 rx3  V  vx 


Ur z 

8  e 


C  uU  ,n 

-  47T  V  2  _  Z 


z  e 


(C.69) 


The  function  (2  —  z)^/ze  z  has  a  maximum  value  of  0.6902,  so  this  term  is  bounded 
by 


E1 


0.69(7  uU 


47T 


x° 


uU 


vU 


—  ~  0.055  CW—  <  0.055  Ca  —r  =  0.025 


pU2r2a 


x 


5  — 


XX 


XX 


(CTO) 
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(n\ 

Since  \vx  |  is  largest  on  the  centreline  (see  Figs.  A. 2,  A. 5  and  A. 8)  and  increases 

(c2'\ 

monotonically  with  x  (see  Figs.  A. 3,  A. 6  and  A. 9),  we  can  estimate  Vx  ;  by  evaluating 
Eq.  (C.44)  with  z  —  0  and  x  =  xmi(i  =  1.5x0.  For  x  <  xmid, 


C2 

v{2)  <  -1.1  x  10~3 - =  -5.5  x  KT3/?2!/ 

VXq 

Syn  pp 

dvf]  CUr  _uS  C  [Uz'  _z 

dx  8nux2  47t  V  vx3 

and  since  \[ze~z  has  a  maximum  of  0.4289, 


(C.71) 


(C.72) 


rh/1) 

9  <  0.4289  — 


U 


<  0.034  Ci 


U  (3U 

—  =  0.076  — 


dx  47T  V  VX3  \l  UXn  Xo 

Therefore  the  second  term  on  the  right  hand  side  of  Eq.  (C.68)  is  of  order 


(C.73) 


E2 


-v. 


A,/1) 

(2)  ove 
dx 


4/?3!/2 


-4.2  x  10 


x0 


(C.74) 


From  Eqs.  (C.10)  and  (C.60),  the  third  term  on  the  right  hand  side  of  Eq.  (C.68)  is 


(i) 


*3  = 


c3 


2567T3  V  v3xb 


U 


m  - 2(1  -«-)«-(*■*,) + £ 


ZX0 


n= 2 


Xh 


X 


\fz  (  2e 


(1  —  e~z) 


(C.75) 


For  x  <  xmid,  the  complicated  term  in  z  has  a  maximum  of  about  0.4  occurring  on 
the  initial  plane.  Therefore  for  x  <  xmid  we  have 


E3  <0.4 


C3 


u 


2567T3  V  ^3x5 


<  5  x  lO'5^3 


JJ  R'i  T  T'2 


3^5 


6  x  10 


4 FU2 


v°x\ 


x0 


(C.76) 


Since  for  the  values  of  (3  used  in  our  calculations,  (3  rCi/x o,  E\  is  significantly 
smaller  than  E2  and  E3,  while  E2  and  E3  are  of  similar  size.  Therefore  we  can 
approximate  the  magnitude  of  the  right  hand  side  of  Eq.  (C.68)  by  E3.  Then  we  can 
estimate  that 


dtj3)  E3 
dx  U 


6  x  10~4 


(33U 


x0 


(C.77) 
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and  therefore  at  x  =  xmid 


(C.78) 


The  maximum  value  of  v j}1  at  the  mid-plane  is  0.093/3E.  Therefore  the  relative  error 
in  the  approximation  of  Vg  by  ev^  is 


(C.79) 


In  practice  this  is  a  conservative  estimate,  both  because  we  have  used  maximum 
values  when  estimating  the  size  of  the  error  terms,  and  because  we  have  neglected 
the  diffusive  term  in  the  estimate  of  Eq.  (C.77). 

When  f3  =  0.0446,  E  ~  2  x  10-5  and  when  f3  =  0.446,  E  ~  2  x  10-3.  Therefore 
for  these  values  of  f3  the  approximate  solution  is  very  good  and  we  should  expect 
the  numerical  solutions  to  converge  to  it.  When  (3  =  4.46,  E  «  0.2  so  that  we 
should  expect  differences  of  the  order  of  20%  between  the  numerical  and  approximate 
solution.  This  is  roughly  consistent  with  the  results  shown  in  Fig.  A. 7. 

Of  the  terms  on  the  right  hand  side  of  Eq.  (C.67),  the  Erst  is  proportional  to  f3 4 
while  the  others  are  proportional  to  [32tE/xq.  Therefore  the  first  term  is  dominant 
for  the  values  of  (3  that  we  have  used.  Comparing  this  term  to  (u^)2/r,  which  is 
the  equivalent  term  for  the  radial  gradient  of  p(2\  suggests  that  the  relative  error  in 
the  approximation  of  p  is  approximately  twice  the  relative  error  in  vg.  Again,  this  is 
roughly  consistent  with  the  results  shown  in  Fig.  A. 7. 
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List  of  symbols 


a 

(3 


v 


Vt 

C 

Cp 

J 

Ne 

Nc 

Nr 

Nx 

P 

( r ,  0 ,  x) 
rc 


R 


so 


sr 

U 

(vx,Vr,Ve) 

V 

K 


Constant  of  proportionality  relating  ut  and  C  for  a  fully-developed 
turbulent  vortex:  vt  ~  aC. 

Non-dimensionalized  circulation:  C/Urci. 

Kinematic  viscosity. 

Turbulent  kinematic  eddy  viscosity. 

The  circulation  of  the  vortex. 

The  pressure  coefficient:  Cp  —  p/ \pU2. 

Advance  coefficient. 

The  number  of  nodes  around  the  circumference  of  the  hexahedral  grid. 

The  number  of  nodes  across  the  vortex  core. 

The  number  of  nodes  radially  in  the  hexahedral  grid. 

The  number  of  nodes  in  the  x  direction  in  the  hexahedral  and  prismatic 
grids. 

Pressure. 

Cylindrical  coordinates. 

The  radius  of  the  vortex  core. 

The  radius  of  the  vortex  core  in  the  upstream  plane. 

The  outer  radius  of  the  grid. 

The  distance  between  radial  nodes  at  r  =  0  in  the  hexahedral  grid. 

The  distance  between  radial  nodes  at  r  —  R  in  the  hexahedral  grid. 

The  magnitude  of  the  free-stream  velocity. 

Components  of  the  velocity  in  a  cylindrical  coordinate  system. 

Speed  of  advance  of  a  propeller. 

Radial  component  of  velocity  past  a  propeller.  The  radial  direction  is 
relative  to  the  propeller  axis. 
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The  x  coordinate  of  the  upstream  end  of  the  grid. 
The  x  coordinate  halfway  down  the  grid. 

The  x  coordinate  of  the  downstream  end  of  the  grid. 


x0 

Xmid, 

X\ 
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